#include <cmath>
#include "openturns/Debye.hxx"

/** 1/(2pi) */
#define M_1_2PI .159154943091895335768883763373

/** Debye function of order n.
    int(t=0..x) t^n dt / [exp(t)-1]
    The underivative for n=0 is log[1-exp(-x)], which is infinite at x=0,
    so the corresponding Debye-function is not defined at n=0.

    Literature:
    Ng et al, Math. Comp. 24 (110) (1970) 405
    Guseinov et al, Intl. J. Thermophys. 28 (4) (2007) 1420
    Engeln et al, Colloid & Polymer Sci. 261 (9) (1983) 736
    Maximon , Proc. R. Soc. A 459 (2039) (2003) 2807

    @param[in] x the argument and upper limit of the integral. x>=0.
    @param[in] n the power in the numerator of the integral, 1<=n<=20 .
    @return the Debye function. Zero if x<=0, and -1 if n is outside the
    parameter range that is implemented.
    @since 2007-10-31 implemented range n=8..10
*/
double debyen(const double x, const int n)
{
  if (x <= 0. )
    return 0. ;
  if ( n < 1 || n > 20)
    return -1. ;

  /* for values up to 4.80 the list of zeta functions
     and the sum up to k < K are huge enough to gain
     numeric stability in the sum */
  if(x >= 3. )
  {
    double sum ;
    static double nzetan[] = { 0., 1.64493406684822643647241516665e+00, 2.40411380631918857079947632302e+00,
                               6.49393940226682914909602217926e+00, 2.48862661234408782319527716750e+01,
                               1.22081167438133896765742151575e+02, 7.26011479714984435324654235892e+02,
                               5.06054987523763947046857360209e+03, 4.04009783987476348853278236554e+04,
                               3.63240911422382626807143525567e+05, 3.63059331160662871299061884284e+06,
                               3.99266229877310867023270732405e+07, 4.79060379889831452426876764501e+08,
                               6.22740219341097176419285340896e+09, 8.71809578301720678451912203103e+10,
                               1.30769435221891382089009990749e+12, 2.09229496794815109066316556880e+13,
                               3.55688785859223715975612396717e+14, 6.40238592281892140073564945334e+15,
                               1.21645216453639396669876696274e+17, 2.43290316850786132173725681824e+18,
                               5.10909543543702856776502748606e+19, 1.12400086178089123060215294900e+21
                             } ;

    /* constrained to the list of nzetan[] given above */
    if ( n >= static_cast<int>(sizeof(nzetan) / sizeof(double)) )
      return -1. ;

    /* n!*zeta(n) is the integral for x=infinity , 27.1.3 */
    sum = nzetan[n] ;

    /* the number of terms needed in the k-sum for x=0,1,2,3...
     * Reflects the n=1 case, because higher n need less terms.
     */
    static int kLim[] = {0, 0, 0, 13, 10, 8, 7, 6, 5, 5, 4, 4, 4, 3} ;

    const int kmax = (static_cast<int>(x) < static_cast<int>( sizeof(kLim) / sizeof(int))) ? kLim[static_cast<int>(x)] : 3 ;
    /* Abramowitz Stegun 27.1.2 */
    for(int k = 1; k <= kmax ; k++)
    {
      /* do not use x(k+1)=xk+x to avoid loss of precision */
      const double xk = x * k ;
      double ksum = 1. / xk ;
      double tmp = n * ksum / xk ;  /* n/(xk)^2 */

      for (int s = 1 ; s <= n ; s++)
      {
        ksum += tmp ;
        tmp *= (n - s) / xk ;
      }
      sum -= exp(-xk) * ksum * pow(x, n + 1.) ;
    }
    return sum ;
  }
  else
  {
    static double koeff[] = { 0., 1.289868133696452872944830333292e+00, 1.646464674222763830320073930823e-01,
                              3.468612396889827942903585958184e-02, 8.154712395888678757370477017305e-03,
                              1.989150255636170674291917800638e-03, 4.921731066160965972759960954793e-04,
                              1.224962701174096585170902102707e-04, 3.056451881730374346514297527344e-05,
                              7.634586529999679712923289243879e-06, 1.907924067745592226304077366899e-06,
                              4.769010054554659800072963735060e-07, 1.192163781025189592248804158716e-07,
                              2.980310965673008246931701326140e-08, 7.450668049576914109638408036805e-09,
                              1.862654864839336365743529470042e-09, 4.656623667353010984002911951881e-10,
                              1.164154417580540177848737197821e-10, 2.910384378208396847185926449064e-11,
                              7.275959094757302380474472711747e-12, 1.818989568052777856506623677390e-12,
                              4.547473691649305030453643155957e-13, 1.136868397525517121855436593505e-13,
                              2.842170965606321353966861428348e-14, 7.105427382674227346596939068119e-15,
                              1.776356842186163180619218277278e-15, 4.440892101596083967998640188409e-16,
                              1.110223024969096248744747318102e-16, 2.775557561945046552567818981300e-17,
                              6.938893904331845249488542992219e-18, 1.734723476023986745668411013469e-18,
                              4.336808689994439570027820336642e-19, 1.084202172491329082183740080878e-19,
                              2.710505431220232916297046799365e-20, 6.776263578041593636171406200902e-21,
                              1.694065894509399669649398521836e-21, 4.235164736272389463688418879636e-22,
                              1.058791184067974064762782460584e-22, 2.646977960169798160618902050189e-23,
                              6.617444900424343177893912768629e-24, 1.654361225106068880734221123349e-24,
                              4.135903062765153408791935838694e-25, 1.033975765691286264082026643327e-25,
                              2.584939414228213340076225223666e-26, 6.462348535570530772269628236053e-27,
                              1.615587133892632406631747637268e-27, 4.038967834731580698317525293132e-28,
                              1.009741958682895139216954234507e-28, 2.524354896707237808750799932127e-29,
                              6.310887241768094478219682436680e-30, 1.577721810442023614704107565240e-30,
                              3.944304526105059031370476640000e-31, 9.860761315262647572437533499000e-32,
                              2.465190328815661892443976898000e-32, 6.162975822039154730370601500000e-33,
                              1.540743955509788682510501190000e-33, 3.851859888774471706184973900000e-34,
                              9.629649721936179265360991000000e-35, 2.407412430484044816328953000000e-35,
                              6.018531076210112040809600000000e-36, 1.504632769052528010200750000000e-36,
                              3.761581922631320025497600000000e-37, 9.403954806578300063715000000000e-38,
                              2.350988701644575015901000000000e-38, 5.877471754111437539470000000000e-39,
                              1.469367938527859384580000000000e-39, 3.673419846319648458500000000000e-40,
                              9.183549615799121117000000000000e-41, 2.295887403949780249000000000000e-41,
                              5.739718509874450320000000000000e-42, 1.434929627468612270000000000000e-42
                            } ;

    double sum = 0. ;

    /* Abramowitz-Stegun 27.1.1 */
    const double x2pi = x * M_1_2PI ;
    for(int k = 1; k < static_cast<int>(sizeof(koeff) / sizeof(double)) - 1 ; k++)
    {
      const double sumold = sum ;
      /* do not precompute x2pi^2 to avoid loss of precision */
      sum += (2. + koeff[k]) * pow(x2pi, 2.*k) / (2 * k + n) ;
      k++ ;
      sum -= (2. + koeff[k]) * pow(x2pi, 2.*k) / (2 * k + n) ;
      if(sum == sumold)
        break ;
    }
    sum += 1. / n - x / (2 * (1 + n)) ;
    return sum * pow(x, (double)n) ;
  }
}

#undef M_1_2PI
